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ABSTRACT: Bolometric detectors are used in particle physics experiments to search for rare pro- 
cesses, such as neutrinoless double beta decay and dark matter interactions. By operating at cryo- 
genic temperatures, they are able to detect particle energies from a few keV up to several MeV, mea- 
suring the temperature rise produced by the energy released. This work focusses on the bolometers 
of the CUORE experiment, which are made of Te02 crystals. The response of these detectors 
is nonlinear with energy and changes with the operating temperature. The noise depends on the 
working conditions and significantly affects the energy resolution and the detection performances 
at low energies. We present a software tool to simulate signal and noise of CUORE-like bolome- 
ters, including effects generated by operating temperature drifts, nonlinearities and pileups. The 
simulations agree well with data. 
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1. Introduction 

Bolometers are detectors in which the energy from particle interactions is converted to heat and 
measured via their rise in temperature. They provide excellent energy resolution, though their 
response is slow compared to electronic or photonic detectors. These features make them a suit- 
able choice for experiments searching for rare processes, such as neutrinoless double beta decay 
(OvDBD) and dark matter (DM) interactions. 

The CUORE experiment will search for OvDBD of 130 Te [jl|, g] using an array of 988 Te0 2 
bolometers of 750 g each. It may also be sensitive to DM interactions [§J. Operated at a temperature 
of about lOmK, these detectors exhibit an energy resolution of a few keV over an energy range 
extending from a few keV up to several MeV. In this range the response function is found to be 
nonlinear The conversion from signal amplitude to energy is complicated and the shape of 
the signal depends on the energy itself. Moreover, the amplitude of the signal depends on the 
temperature of the detector, which is very difficult to keep stable with current cryostats within the 
few ppm level, the level that would not perturb the energy resolution. The noise of the detector 
is dominated by thermal fluctuations induced by vibrations, and significantly affects the energy 
resolution at low energies [Qj. 

In this paper we present a method to simulate signal and noise of CUORE-like bolometers. The 
simulation should be able to reproduce all the features of the data and can be used, for example, 
to estimate detection efficiencies and to test analysis algorithms. The shape of the signal and the 
noise are estimated from the data. The nonlinearities of the signal are reproduced using a model of 
the thermal sensor of the bolometer [fflj. 
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2. Description of the detector 



A CUORE bolometer is composed of two main parts, a Te02 crystal and a neutron transmutation 
doped Germanium (NTD-Ge) thermistor []6[ ^]. The crystal is cube-shaped (5x5x5 cm 3 ) and held 
by Teflon supports in copper frames. The frames are coupled to the mixing chamber of a dilu- 
tion refrigerator, which keeps the system at a temperature of ~ 10 mK. The thermistor is glued to 
the crystal and acts as thermometer (Fig. [[]). A Joule heater is also glued to most crystals. It is 
used to inject controlled amounts of energy into the crystal, to emulate signals produced by parti- 
cles [M, 9fl. When energy is released in the crystal, the crystal temperature increases and changes 



the thermistor's resistance according to the relationship [10]: 



R(T)=R exp(T /Ty 



(2.1) 



where Ro and 7b are parameters that depend on the dimensions and on the material of the thermistor. 
For CUORE bolometers values are about 1.1 O and 3.4 K, respectively. At 10 mK the parameter y 
can be considered constant and equal to 1/2 [jlT, 12]. 



////////// • 
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Figure 1. Sketch of a CUORE-like bolometer (left) and a photograph of a bolometer (right). The TeC>2 
crystal is held by Teflon supports, the thermistor is glued to the crystal and its wires are attached to the 
copper frame. The supports and the thermistor wires thermally couple the crystal to the copper frame, which 
act as heat bath. 



To read out the signal, the thermistor is biased in differential configuration with a bipolar volt- 
age generator ±Vbias connected to a pair of load resistors, Rl's, finally connected to the thermistor's 
terminals. The resistance of the thermistor varies in time with the temperature, R(t), and the voltage 
across it, V#(f), is the bolometer signal. The value of R^s is chosen to be much higher than R(t) 
so that Vr(?) is proportional to R(t). Since from Eq. (2.1) positive temperature variations induce 
negative resistance variations, the polarity of Vbias is chosen to be negative in order to obtain pos- 
itive signals. The connecting wires add in parallel to the thermistor a parasitic capacitance c p . A 
schematic of the biasing circuit is shown in Fig. ||. In the figure, the series of the two load resistors 
is represented as a unique resistor R^. The signal Vs(?) is amplified, filtered with a 6-pole active 
Bessel filter, and then digitized with an 18-bit analog-to-digital converter (ADC). To fit the signal 
in the range of the ADC, which is [-10.5,10.5] V, a programmable offset voltage V/, is added to 
Vjt(t). The front-end electronics, which provide the bias voltage, the load resistors, the amplifier 
and the offset voltage, are placed outside of the cryostat, at ambient temperature [13]. 

At lOmK the value of R(T) is of order 100MQ, R L is chosen as 54GQ. (27 + 27 GO.) and 
Vbias as ~ 5 V. The c p value depends on the length of the wires that carry the signal out of the 
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Rl 




Figure 2. Biasing circuit of the thermistor. A voltage generator Vbias biases the thermistor resistance R(t) 
in series with a load resistance Rl. The bolometer signal is the voltage Vr(?) across R(t ). The wires used to 
extract Vg(?) from the cryostat have a non-negligible capacitance c p . 



cryostat, typically it is of order 400pF. The amplifier gain, the Bessel filter frequency bandwidth, 
the duration of the acquisition window and the sampling frequency are set typically at 5000 V/V, 
12Hz, 5.008s and 125 Hz, respectively. 

The data analyzed in this paper come from test bolometers operated by the CUORE collabora- 
tion at the Gran Sasso underground laboratory (LNGS) in Italy [|l4|]. The bolometers were exposed 
to a 232 Th calibration source which, together with an a line generated by 210 Po contamination in 
the crystal, allows the analysis of an energy range up to 5407 keV. 

To simplify the description of this work we will focus our analysis on a single bolometer. The 
signal rate on that bolometer was 133 mHz, that has to be combined to the rate of heater pulses that 
were fired at an energy of 1885keV every 300 seconds (3.3 mHz). The energy spectrum acquired 
in about 3 days is shown in Fig. ||. 




Energy (keV) 

Figure 3. Energy spectrum. All lines are generated by the 232 Th calibration source except for the line at 
5407 keV, arising from 210 Po contamination in the TeC>2 crystal. Heater pulses were fired at an energy of 
1885keV. 
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3. Signal and noise features 



Examples of signals generated by a 2615keV /-ray and by the heater, as acquired by the ADC, are 
shown in Fig. ||. The baseline voltage of the pulses is related to the thermistor temperature in static 
conditions, and the amplitude is related to the energy released. The shape of heater pulses is found 
to be slightly different from that of particle pulses: the rise time of the particle (heater) pulse in the 
figure, computed as the time difference between the 10% and the 90% of the leading edge, is 55 
(54) ms while the decay time, computed as the difference between the 90% and 30% of the trailing 
edge, is 220 (255) ms. 
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Figure 4. Pulse shapes of a 2615 keV y-ray (left) and heater (right). The baseline is related to the temperature 
of the thermistor before the particle interaction or the heater shot. The amplitude carries information on the 
amount of energy released. 



As already observed in Ref. [Q] several nonlinearities are present: 
1. The rise and the decay times of a pulse depend on the energy (Fig. ||). 



2. The amplitude of the pulse depends on the base temperature, which varies during the data 
acquisition (Fig. ||). 




Amplitude (mV) Amplitude (mV) 

Figure 5. Pulse shape parameters versus energy. The correlation with energy is negative for the rise time 
(left) and positive for the decay time (right). Heater pulses are marked in red. 
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Figure 6. Amplitude of heater pulses versus baseline. A change in the bolometer temperature also changes 
its response, degrading the energy resolution. 
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Figure 7. Residuals obtained using a linear calibration function for the well-identified peaks in the 232 Th 
source spectrum. Considering that the energy resolution is ~ 5keVFWHM, the difference AE between 
the estimated energies and the true peak energies is not compatible with zero. The error bars refer to the 
uncertainty on the estimated peak position that depends on the FWHM and on the number of events in the 
peak as FWHM/(2.35V^V). 



3. The amplitude dependence on energy is not linear. The deviation from linearity of the data 
is estimated by comparing the result of a linear calibration function, 

Energy = constant • Amplitude , (3.1) 

to the true energy of the source peaks. The residuals evaluated on the peaks generated by 
the 232 Th source (see Fig. [3]) are shown in Fig. [7|. The 5407 keV line is not considered here 
because a particles have a quenching factor different from / and j8 particles [Q]. 

The noise of the bolometer in the signal frequency region (0-10 Hz) is a result of vibrations 
inside the cryostat, and depends on the mechanical setup of the experimental apparatus. The noise 
power spectrum was estimated as: 

N((O k ) =<\n(co k )\ 2 > (3.2) 
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where n((Ok) is the k-th component of the discrete Fourier transform (DFT) of an acquired wave- 
form not containing signals, and <> denotes the average over a large number of waveforms. A 
typical noise waveform and the estimated power spectrum are shown in Fig. ||. The peaks in the 
power spectrum are due to the crystal friction against its frame, and the residual common mode 
contribution from the vibration of the connecting wires (readout by the differential preamplifier). 
Which of the two sources is dominant may depend on the set-up. The white noise contribution at 
high frequency is due to the ADC digitization. 

In the next sections we will describe the procedure we developed to simulate the signal shape, 
the nonlinearities and the noise. 




Time (s) Frequency [Hz] 

Figure 8. Noise of a bolometer at the output of the acquisition chain. A sample waveform (left) and the 
power spectrum estimated from a large number of waveforms (right). 



4. Signal model 

To simulate the data we developed a model of the bolometer signal that reproduces the pulse shape 
and the nonlinearities. Starting from the energy release in the crystal, the model is divided into 
three main stages: the thermal response of the bolometer, the response of the thermistor and of its 
biasing circuit, and the simulation of the electronics. All parameters are detector variables, except 
for the thermal model parameters, which we determine from fits to the data. 

4.1 Thermal model 

A bolometer is a thermal system composed of a crystal, crystal supports, thermistor, and their cou- 



pling elements [ ]15[ ] (Fig. |9|). The crystal capacitance C c is connected to the thermistor through 
the glue spots with conductance K g and to the supports through a contact conductance K cs . The 
thermistor can be represented as a two-stage system composed of a lattice and an electron gas, each 
with its capacitance and conductance, and inter-connected by a conductance K ep . The lattice ca- 
pacitance, not shown in the figure, is negligible and it discharges through the gold wires connected 
to the main heat bath Ka u . The electron gas capacitance and conductance are labeled as C e and K e , 
respectively. The left side of the circuit shows the crystal supports with their capacitance C s and 
their conductance to the main bath K s . The main bath acts as the reference ground. 

We are interested in the expression of the temperature variation of the thermistor electron gas 
(node 4 in Fig. |9|) as a function of the time AT(t) after that an amount of energy E is released in 
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Figure 9. Thermal circuit of a CUORE-like bolometer. Each of the elements is identified in the text. 



the crystal. The energy release is effectively instantaneous, because the phonons produced by a 
particle thermalize in a time much shorter than the rise time of acquired pulses JTr^]. Under these 
conditions the analytical expression of AT(t) is found to be: 

AT(t) =A T (-e~* +ae - vT + (l-a)e~^ ; (4.1) 

where all the parameters are complicated functions of the elements of the thermal circuit. The 
formula indicates that the thermistor temperature increases with one rise time constant, x r , and 
decreases with two decay constants, T^i and Tj2- The parameter a weighs the two exponential 
decays and satisfies the condition < a < 1. The amplitude of the thermal pulse, A 7-, is also a 
function of the thermal elements and is directly proportional to E. 

In principle the parameters in Eq. [O could be known if one were able to measure the un- 



derlying thermal elements. Previous investigators have measured the thermal parameters [|15[], and 
obtained heat capacitances and heat conductances of order 10~ a J/K and 10~ fo W/K with a in the 
range 9 to 10 and b in the range 9 to 11. Such measurements are of little utility for us for several 
reasons. First, the capacitances and conductances depend on the temperature and should be mea- 
sured in the bolometer working temperatures, that depend on the setup. Second, several thermal 
elements vary from bolometer to bolometer. For example, the contact heat conductance between 
the crystal and its supports (K cs ) changes with the detector configuration, and the heat conduc- 
tance of the glue (K g ) varies because of the weak reproducibility of the glue deposition. Third, we 
should know the parameters with a precision of the order of the present energy resolution (0.1%), 
which is not feasible using ordinary measurements techniques. We choose instead to estimate the 
parameters in Eq. 4J from fits to the signal waveforms (see Sec. |5]). 



The solution of the circuit as expressed in Eq. 4.1 does not include the dependence of the 
thermal elements on the temperature. This potential source of nonlinearity, however, is not visible 
in the data [ffl]. We also neglected the electrothermal feedback effect: when temperature change 
induces resistance change in a biased thermistor, the resulting excursion in Joule heating induces 
further temperature change. It can be shown that the electrothermal feedback to the first order acts 
as a correction to the values of K e and C e , and therefore does not affect the configuration of the 



thermal circuit and the form of Eq. 4. 1 . 
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4.2 Thermistor and biasing circuit models 



When the thermistor temperature varies, its resistance varies according to Eq. [2T 

T 



AR(AT) =R exp 



T B + AT 



R* 



(4.2) 



where T B is the initial temperature of the bolometer and R B = R(T B ). Because we do not know the 
parameters Rq, Tq and T B with the desired accuracy, we adopt an approximation valid for the small 
AT that will occur [Qj: 

AR(AT) ~/? B [exp(-T]Ar/r B )-l] , (4.3) 

where 



dlogR 



dlogT 



7 log 



R(T) 



Rn 



(4.4) 



v\ is the sensitivity of the thermistor, which has value of order 10 but is not known with precision. 
The advantage of Eq. 13 is that the parameter R B can be measured with precision and that the un- 
known f]/T B is just a scale factor applicable to the temperature variation. We obtain the expression 



for the resistance variation after an energy release by substituting Eq. [kjj into Eq. ^3 
AR(t) = R B jexp 



-e +ae vi + (l - a)e ^ 



1 



(4.5) 



where A = T}At/T b . This expression absorbs the unknown parameter f]/T B with the unknown 
thermal amplitude A j. The thermistor model therefore does not change the number of unknown 
parameters and adds the measurable parameter R B . 

The relationship between the voltage across the thermistor V R {t) and its resistance R(t) can be 
obtained from the differential equation describing the thermistor's biasing circuit (see Fig. ||): 



R L +R(t) 
R(t) 



V R (t)-V b ias+RLCp 



dV R (t) 
dt 



0. 



(4.6) 



The model we are building is based on variations of the resistance from the measured value of R B , 
which in turn generate voltage variations from the corresponding voltage V B : 



V R = Vbias- 



R b 



'R B + R L ' 

By splitting R(t) and Vg(t) into time-independent and time-dependent contributions, 

R(t)=R B + AR(t) V R (t)=Vg + AV R (t), 
we obtain the differential equation relating resistance and voltage variations: 



R L + R B + AR(t) 
R B + AR(t) 



Vbias 



R b 



R b + Rl 



+ AV R (t) 



dAV R (t) 
V bias +R L c p — —=0 
dt 



(4.7) 



(4.8) 



(4.9) 



Given the form of AR(t) in Eq. [4.5| , AV R (t) cannot be obtained in an closed form. In our bolometer 
model we solve Eq. 4.9 numerically using the Runge-Kutta method [|l7|]. 
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The thermistor and the biasing circuit are the only sources of nonlinearities of our model, and 
should be able to describe the nonlinearities observed in the data. If we assume that the thermal 
circuit responds linearly, the amplitude A in Eq. is directly proportional to the energy E released 
in the bolometer, 

A = c-E. (4.10) 
The corresponding resistance variation, however, is not proportional to A, because the exponential 



dependency in Eq. |45| does not transform linearly the shape and the amplitude of the pulse. More- 
over the voltage variation is not strictly proportional to the resistance variation. This model descrip- 
tion should be sufficient to generate the shape dependence on energy and the nonlinear calibration 
function (Figs. || and [7]). The amplitude dependence on the baseline in Fig. ^ can be generated from 
Eq. by varying the baseline voltage of the pulse and hence the thermistor resistance R B . 

4.3 Electronics 

The front-end electronics amplifies the bolometer signal AVg(t) by G, a parameter that is measured 
with a precision better than 0.2%. The deviation from linearity of the amplifier, in the voltage range 
we are considering, is less than 0.01%. The output voltage of the amplifier 

AV G (t)=AV R (t)-G, (4.11) 

is fed into a six-pole Bessel filter, whose transfer function is 

10395 

^ ~ a 6 + 21a 5 + 210a 4 + 1260a 3 +4725a 2 + 10395a + 10395 ' 

In the above equation a is the normalized Laplace variable which can be expressed in terms of the 
frequency ft) as: 

2.703395061 

d=Jt0 , (4.13) 

Jb 

where ft is the filter cutoff (12 Hz in our case). 

The signal is filtered by multiplying its DFT, AVg(co), by fi(ft)), removing the "DFT wraparound 
problem" with the method described in Ref. JT8|]. The output of the filter is then obtained as: 

AV(t) = &- l [AV G {(D) -fi(ft))] , (4.14) 

where J£"~ 1 denotes the inverse DFT. 



In summary , the model of the signal, from the energy release in the crystal E to the signal 
acquired by the ADC AV(t), is obtained using the thermal model in Eq. the thermistor model 
in Eq. [O], the voltage across the thermistor from Eq. [4^ , the amplifier and Bessel filter effects in 
Eqns. |TT] and |414| : 



AV(t) 



Eq:M 
analytic 
solution 



>AT(t) 



Eq. 
analytic 
solution 



^AR(t) 



numerical 
solution 



*AV R (t) 



Eq. 4.11 



G 



Eq 



B(t). (4.15) 



convolution 
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The baseline voltage of the signal, V B , is the sum of two components, the thermistor voltage 
V B in Eq. 4/7 scaled by the electronics gain G, and the offset voltage Vf, added by the electronics 
itself: 



V B = v£-G + V h = V bias G 



R b +Rl 



+v h . 



(4.16) 



To reproduce a real waveform, V(t), we add the baseline voltage to the pulse in Eq. 4.15 and we 
also account for the onset time to of the pulse, which starts about 1 s after the beginning of the 
waveform: 



V(t) = V B + &(t-t )AV(t-t ) 



(4.17) 



where &(t) is the Heaviside step function. Since the thermistor temperature is not stable, the 
resistance R B and the baseline voltage V B are not fixed parameters of the model. V B is measured on 
data by averaging the first 0.8 s of the waveform. The corresponding value of R B is then computed 
from Eq. 4.16| and used in the signal model. In Tab. [j] we list the parameters of the model and 
indicate whether they are measured, or determined from fits to signal waveforms. 



Table 1. Parameters of the signal model in Eqns. |4.1 
Parameter Name 



4.16 and 4.17 



Equation Estimation 



It 


Thermal rise time 


4.1 


Fit 


a 


Weight of the two thermal decay constants 


4.1 


Fit 




Fast thermal decay constant 


4.1 


Fit 


Idl 


Slow thermal decay constant 


4.1 


Fit 


c 


Energy to thermal amplitude conversion 


4.10| Fit 


R B 


Thermistor resistance at the pulse baseline 


4.5 


Measured 


Vbias 


Bias voltage 


4.9 


Measured 


Rl 


Load resistor 


4.9 


Measured 


c p 


Parasitic capacitance 


4.9 


Measured 


G 


Electronics gain 


4.11 


Measured 


fb 


Bessel filter cutoff frequency 


4.U 


Measured 


v h 


Electronics offset voltage 


4.16 Measured 


V B 


Baseline voltage 


4.17 Measured 


to 


Onset time of the pulse 


4.17| Fit 



5. Estimation of the signal model 

We intend that the model we developed accounts for all the nonlinearities of the signal. The un- 
known thermal parameters in Tab. [j] are expected to be independent of the energy, except for the 
amplitude A, which should be proportional to the energy (Eq. 4.10j ). The energy independence 



allows us to determine unmeasured parameters from fits to pulses at a single energy, and apply the 
resulting model over the entire range of energies. We performed fits on particle pulses occurring in 
an energy window of 30keV around the 2615keV 7 peak. We performed a separate set of fits on 
heater pulses, because their shape differs from the shape of particle pulses. 
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Figure 10. Fit of a 2615keV pulse. Data (black solid lines) with superimposed fit function (Eq. 4.17) (blue 
dashed lines) are shown in the left column, fit residuals are shown in the right column. The two bottom rows 
display a zoom of the rise and of the maximum of the pulse, respectively. The fit region ends when the pulse 
falls to within 3 standard deviations of the baseline level. 



Table 2. Average parameters fitted on particle (2615keV 7) and heater pulses and the j 2 /ndf of the fits. 



Parameter 


Particle 


Heater 


x r (ms) 


20.70 ±0.05 


18.8 ±0.1 


a 


0.916 ±0.002 


0.911 ±0.002 


x dX (ms) 


158.1 ±0.5 


187.3 ±0.6 


T d2 (ms) 


770 ±20 


970 ±20 


c (1/MeV) 


0.04703 ±0.00006 


0.04356 ±0.00009 


to (s) 


1.0145 ±0.0005 


0.9944 ±0.0005 


Z 2 /ndf 


2.4 ±0.3 


2.0 ±0.2 



Figure [10| shows a typical fit to a 2615keVy-ray pulse, Tab. ^reports the parameters averaged 
over 25 fits of particle and heater pulses. The average % 2 is about twice the number of degrees 
of freedom, i.e. about twice the value one expects for a perfect model. Shortcomings of the model 
are also apparent in the fit residuals on the right of the figure, where mismatches are evident in 
the leading edge and in the vicinity of the maximum of the pulse. The voltage amplitude of the fit 
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function is slightly biased, and found to be, on average, higher than the amplitude of the pulse by 
0. 17 ± 0.02 % for particle pulses and 0.37 ± 0.02 % for particle heater pulses. 

The fit error can be ascribed to an incompleteness of the thermal model, that would probably 
benefit from the inclusion of second order effects like the nonlinearities of the thermal elements 
and the electrothermal feedback. For our purposes, however, the model performs well, reproducing 
the signal at the per mil level. 



6. Noise generation 

A complete and predictive model for the noise of CUORE bolometers is missing. The main sources 
of fluctuations that spoil the energy resolution have different origins and, to estimate the overall 
contribution, each of them should be propagated with the transfer function of each step of the 
acquisition chain. Examples of these sources are the Johnson noise of the load resistors, vibrations 
of the experimental apparatus that dissipate energy in the bolometer and instabilities of the cryostat 
temperature. All these effects contribute to the power spectrum shown in Fig. || 

Although the noise from the load resistors and amplifiers is predictable, a model for the noise 
from vibrations is not in hand. Moreover, it varies from bolometer to bolometer, because of the 
weak reproducibility of the assembly. On this account we adopted a statistical model for the noise. 
We simply require that, on average, the simulated time series behave like the experimental one, 
namely that the average power spectrum of the simulated baselines is as close as possible to data. 



We followed an approach [ |19| ] based on an application of the Carson's theorem Q20| ] to a 
discrete time series. A random waveform n(//) can be represented as a superposition of independent 
pulses of fixed shape g(t{) and amplitude A, distributed in time according to a Poisson process of 
rate X : 

n{ti) =AY d g(t i -ti) , (6.1) 

/ 

where the differences between consecutive f/'s follow an exponential distribution with mean 1/X. 
The theorem states that N(cOjc) and G(cOk), the average power spectra of n(tj) and respectively 



(see Eq. |3.2[ ), satisfy the relationship: 

N(co k ) = XTA 2 G(co k ) , (6.2) 

where T is the length of the time series (5.008 s in our case). To utilize the theorem one must find 
a shape g(tj) for which the power spectrum is proportional to the average power spectrum of the 



noise time series to be simulated. A 2 and X are adjustable parameters with product fixed by Eq. ^2. 
These parameters set the aspect of the noise in the time domain, i.e. the same power spectrum can 
be produced with a small rate of pulses with large amplitude and vice versa. Once g(tj) and, say, X 



are determined the method is fully specified and «(?,) is obtained from Eq. |6T 
We built g(tj) as the inverse DFT of 



g(c0k) = ^/N{(Q k y e * , (6.3) 

where d k is a phase randomly sampled within [0, 2n] . The reality of g(tj) is guaranteed by imposing 
the constraint g((Ok) = g*(—Q)k)- X was chosen equal to the Nyquist frequency (62.5 Hz), i.e. the 
maximum rate producing distinguishable pulses. 
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The left panel of Fig. |T1| shows a simulated and a measured noise waveform, and their appeare- 
ance is similar. The right panel of the figure shows the agreement between simulated and measured 
power spectra. 




Time (s) Frequency [Hz] 

Figure 11. Comparison of noise simulation (solid red line) with data (dashed black line / dots). Sample 
waveforms (left) and the power spectrum estimated from a large number of waveforms (right). 



7. Simulation and validation with data 

We built a simulation engine that is able to generate particle, heater, and pure noise waveforms. It 
uses as input the parameters of the signal model in Tab. [I] and the measured noise power spectrum. 
To correct the small error of the model in reproducing the pulse shape (see Sec. |J), the engine scales 
the simulated signals by the estimated bias, so as to match the amplitude of measured signals. The 
simulated noise is summed to the signal assuming that it is purely additive. The additivity is 
supported by the fact that the absolute energy resolution for heater pulses is independent of the 
energy release, and by the residuals in Fig. that do not seem correlated with the time evolution 
of the pulse. 

The simulation can be highly customized. One can choose the energy distribution of the events 
to be generated, the baseline distribution, and the distribution of the time interval between events. 
The features of the signal described in Sec. f| should be automatically reproduced by the model, 
and signals close in time are summed to reproduce pileups. 

To validate the simulation engine, we show the results of a simulation configured to reproduce 
the data shown in Sec. |[ Signals were generated sampling their energy from the spectrum in 
Fig. |5[ the baseline was generated in the range of Fig. |6[ the time delay between particle pulses 
was generated following an exponential distribution with mean 1/(133 mHz), and heater pulses 
were generated every 300 s. Particle and heater pulses were generated using the fitted parameters 
in Tab. ||, and the noise was generated from the power spectrum in Fig. |8[ 

The comparison of the simulation with the data shows good agreement. The shape of the 



pulses with noise added reproduces well the acquired waveforms (Fig. 12). The distributions of the 



rise and decay times shown in Fig. |13] also confirm the agreement. We attribute a small mismatch 
in the rise time of heater pulses to imperfections of the model (see Sec. |5]). The correlation between 
pulse amplitude and baseline is very well reproduced (Fig. Il4|). The dependence of amplitude 



on energy follows the data well at high energies (Fig. 15). The mismatch of order 1% at zero 
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energy is another manifestation of the imperfections of the model. In most applications this error 
is ignorable, nonetheless we incorporated an option to generate waveforms by sampling from the 
amplitude spectrum instead of the energy spectrum, in which case the mismatch is suppressed by 
fiat. 




1 2 3 4 5 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 



Time (s) Time (s) 

Figure 12. Comparison of the simulation (solid red line) of a 2615keV /pulse with data (black dots). Full 
waveform (left) and expansion of the leading edge of the signal (right). 




Amplitude (mV) Amplitude (mV) 



Figure 13. Comparison of the rise time (left) and decay time (right). The black, gray, red, and blue dots 
represent measured particle pulses, simulated particle pulses, measured heater pulses, and simulated heater 
pulses respectively. 
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Figure 14. Comparison of the amplitude of heater pulses versus baseline from data (solid black) and simu- 
lation (open gray). 
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Figure 15. Difference between the amplitudes of simulated and real pulses as a function of the energy (left), 
and deviation from linearity of the calibration function (right) for data (black circles) and simulation (red 
triangles). 



-15- 



8. Applications 



The simulation can be used to test and tune analysis algorithms, comparing the results with the so 
called "Monte Carlo truth", i.e. the generated values of baseline, amplitude, energy and the time 
of each signal. In same cases, in fact, the signal identification is not straightforward and the energy 
could be wrongly estimated. 

For example the data analysis can be complicated by pileups, that alter the baseline and the 
shape of the signals. A simulated sequence of close signals is shown in Fig. 16, where it can be 
seen how a signal can be modified by other signals, and how its identification is complicated. In 
this case the simulation can be used to improve the analysis algorithms and to estimate the error on 
the results. 

A potential application is the estimation of detection efficiencies at low energy. Recently it has 
been demonstrated [^] that CUORE may have an energy threshold of only a few keV, thus being 
sensitive to Dark Matter interactions and rare nuclear decays. Figure [n| shows simulations of very 
small pulses. In this regime noise can both mask signal and mimic signal. Simulations produced 
by our engine may be used to test the immunity of analysis algorithms to both types of error. When 
a heater is available it is used to produce controlled pulses and estimate the detection efficiencies. 
The simulation can correct for the difference in pulse shape between particle and heater generated 
energy deposit. With no heater the simulation could serve on its own to estimate the detection 
efficiencies. 
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Figure 16. Simulated pileups. The pulse at ~ 1 s was generated with an energy of 63keV, lies on the tail of 
a previous pulse of 1899keV and is succeeded by a pulse of 263 keV. 
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Figure 17. 12keV (left) and 6keV (right) simulated pulses. 
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